module mod_chem_phases

   use iso_fortran_env, only: dp => real64

   use mod_chem_elements
   use mod_chem_species
   use mod_chem_reactions
   use mod_chem_transport
   use mod_chem_unit
   use stdlib_optval

#include <chem.f90>

   private
   public :: type_chem_phases

   type :: type_chem_phases

      logical :: flg_Tv = .false.
      real(dp) :: T, P, D, Tv
      real(dp), allocatable :: X(:), Y(:)

      type(type_chem_unit) :: units
      character(len=:), allocatable :: name
      integer :: n_elements
      integer :: n_species
      integer :: n_reactions
      real(dp), allocatable :: elements_molar_mass(:)
      real(dp) :: mean_molecular_weight
      type(type_chem_element), allocatable :: elements(:)
      type(type_chem_specy), allocatable :: species(:)
      type(type_chem_reaction), allocatable :: reactions(:)
      type(type_chem_gas_transport) :: gas_transport
   contains

      procedure :: set_state => SetState

      procedure :: get_temperature => GetTemperature
      procedure :: get_pressure => GetPressure
      procedure :: get_density => GetDensity
      procedure :: get_X => GetX
      procedure :: get_Y => GetY
      procedure :: get_molar_density => GetMolarDensity

      procedure :: get_creation_production_rates => GetCreationProductionRates
      procedure :: get_destruction_production_rates => GetDestructionProductionRates
      procedure :: get_net_production_rates => GetNetProductionRates
      procedure :: get_net_production_rates_dx => GetNetProductionRatesDX
      procedure :: get_fwd_rate_progress => GetFwdRateProgress
      procedure :: get_rev_rate_progress => GetRevRateProgress
      procedure :: get_net_rate_progress => GetNetRateProgress

      procedure :: get_mean_molecular_weight => GetMeanMolecularWeight
      procedure :: get_molecular_weights
      procedure :: get_partial_cp_R => GetPartialCpR
      procedure :: get_partial_cp_mole => GetPartialCpMole
      procedure :: get_partial_cp_mass => GetPartialCpMass
      procedure :: get_partial_enthalpies_RT => GetPartialEnthalpies_RT
      procedure :: get_partial_enthalpies_mole => GetPartialEnthalpiesMole
      procedure :: get_partial_enthalpies_mass => GetPartialEnthalpiesMass
      procedure :: get_partial_internal_energies_RT => GetPartialInterEnergies_RT
      procedure :: get_partial_internal_energies_mole => GetPartialInterEnergiesMole
      procedure :: get_partial_internal_energies_mass => GetPartialInterEnergiesMass
      procedure :: get_partial_entropies_R => GetPartialEntropies_R
      procedure :: get_partial_gibbs_RT => GetPartialGibbs_RT

      procedure :: get_partial_viscosity => GetPartialViscosity
      procedure :: get_partial_thermal_conductivity => GetPartialThermalConductivity
      procedure :: get_viscosity => GetViscosity
      procedure :: get_thermal_conductivity => GetThermalConductivity

      procedure :: get_cp_mole => GetCpMole
      procedure :: get_cp_mass => GetCpMass
      procedure :: get_cv_mole => GetCvMole
      procedure :: get_cv_mass => GetCvMass
      procedure :: get_enthalpy_mole => GetEnthalpyMole
      procedure :: get_enthalpy_mass => GetEnthalpyMass
      procedure :: get_internal_energy_mole => GetInterEnergyMole
      procedure :: get_internal_energy_mass => GetInterEnergyMass
      procedure :: get_gamma => GetGamma
      procedure :: get_sound_speed => GetSoundSpeed

   end type type_chem_phases

contains

   subroutine setState(this, T, X, Y, P, H, U, D)
      class(type_chem_phases), intent(inout) :: this
      real(dp), intent(in), optional :: T
      real(dp), intent(in), optional :: X(:)
      real(dp), intent(in), optional :: Y(:)
      real(dp), intent(in), optional :: P
      real(dp), intent(in), optional :: H
      real(dp), intent(in), optional :: U
      real(dp), intent(in), optional :: D

      integer :: l
      logical :: present_temperature_info

      if (present(T)) call SetTemperature(this, T)
      if (present(P)) call SetPressure(this, P)
      if (present(D)) call SetDensity(this, D)

      if (present(X)) call SetX(this, X)
      if (present(Y)) call SetY(this, Y)

      if (present(X) .eqv. present(Y)) then
         print *, "Error: cannot set both X and Y or not set either"
         error stop
      end if

      if (present(X)) this%Y = X_to_Y(this, X)
      if (present(Y)) this%X = Y_to_X(this, Y)

      present_temperature_info = present(T) .or. present(H) .or. present(U)
      !! check the input variables.
      if ((present_temperature_info .neqv. present(P) .neqv. present(D))) then
         print *, "Error: must set either T/H/U, P, D"
         error stop
      end if

      if (present(H)) then
         this%T = ComputeTemperatureFromH(this, H)
      elseif (present(U)) then
         this%T = ComputeTemperatureFromU(this, U)
      end if

      if (present_temperature_info .and. present(P)) then
         this%D = ComputeDensity(this)
      end if
      if (present_temperature_info .and. present(D)) then
         this%P = ComputePressure(this)
      end if

      if (this%flg_Tv) then
      else
         this%Tv = this%T
      end if

      do l = 1, this%n_species
         call this%species(l)%initial()
      end do

   end subroutine SetState

   subroutine SetTemperature(this, T)
      class(type_chem_phases), intent(inout) :: this
      real(dp), intent(in) :: T

      this%T = T

   end subroutine SetTemperature

   subroutine SetPressure(this, P)
      class(type_chem_phases), intent(inout) :: this
      real(dp), intent(in) :: P

      this%P = P

   end subroutine SetPressure

   subroutine SetDensity(this, D)
      class(type_chem_phases), intent(inout) :: this
      real(dp), intent(in) :: D

      this%D = D

   end subroutine SetDensity

   subroutine SetX(this, X)
      class(type_chem_phases), intent(inout) :: this
      real(dp), intent(in) :: X(:)

      this%X = X

   end subroutine SetX

   subroutine SetY(this, Y)

      class(type_chem_phases), intent(inout) :: this
      real(dp), intent(in) :: Y(:)

      integer :: i

      this%Y = Y

   end subroutine SetY

   function X_to_Y(this, X) result(Y)
      class(type_chem_phases), intent(inout) :: this
      real(dp), intent(in) :: X(:)
      real(dp), allocatable :: Y(:)

      this%mean_molecular_weight = sum(X*this%species%molar_mass)
      allocate (Y(this%n_species))
      Y = X/this%mean_molecular_weight*this%species%molar_mass
   end function X_to_Y

   function Y_to_X(this, Y) result(X)
      class(type_chem_phases), intent(inout) :: this
      real(dp), intent(in) :: Y(:)
      real(dp), allocatable :: X(:)

      this%mean_molecular_weight = 1.0_dp/sum(Y/this%species%molar_mass)
      allocate (X(this%n_species))
      X = Y/this%species%molar_mass*this%mean_molecular_weight

   end function Y_to_X

   function GetTemperature(this) result(T)
      class(type_chem_phases), intent(in) :: this
      real(dp) :: T

      T = this%T

   end function GetTemperature

   function GetPressure(this) result(p)
      class(type_chem_phases), intent(in) :: this
      real(dp) :: p

      p = this%P

   end function GetPressure

   function ComputeTemperatureFromH(this, H) result(T)
      class(type_chem_phases), intent(inout) :: this
      real(dp), intent(in) :: H
      real(dp) :: T

      real(dp) :: dh, T0
      real(dp) :: eps = sqrt(epsilon(1.0_dp))

      T0 = 0.0_dp
      T = H/(3.5_dp*R0/this%mean_molecular_weight) !! initial value
      do while (abs(T - T0) >= sqrt(epsilon(1.0)))
         T0 = T
         dh = (this%get_enthalpy_mass(T0 + eps) - this%get_enthalpy_mass(T0 - eps))/(2.0_dp*eps)
         ! print*, T, dh
         T = T0 - (this%get_enthalpy_mass(T0) - H)/dh
      end do

   end function ComputeTemperatureFromH

   function ComputeTemperatureFromU(this, U) result(T) !! K
      class(type_chem_phases), intent(inout) :: this
      real(dp), intent(in) :: U
      real(dp) :: T, TT

      real(dp) :: du, T0
      real(dp) :: eps = sqrt(epsilon(1.0_dp))
      integer :: iteriation

      T0 = 0.0_dp
      T = U/(2.5_dp*R0/this%mean_molecular_weight) !! initial value
      iteriation = 0
      do while (abs((T - T0)) >= sqrt(epsilon(1.0)))
         iteriation = iteriation + 1
         T0 = T
         du = (this%get_internal_energy_mass(T0 + eps) - this%get_internal_energy_mass(T0 - eps))/(2.0_dp*eps)
         ! print*, T, du
         T = T0 - (this%get_internal_energy_mass(T0) - U)/du
         if (iteriation > 50) then
            ! print*, T, T0, du, this%get_internal_energy_mass(T0), U
            print*, "over maximum iterations!!!!!!!!!!!!"
            print*, iteriation, T, T0, du, this%get_internal_energy_mass(T0), U
            TT = U/(2.5_dp*R0/this%mean_molecular_weight)
            du = (this%get_internal_energy_mass(TT + eps) - this%get_internal_energy_mass(TT - eps))/(2.0_dp*eps)
            print*, TT, du, (this%get_internal_energy_mass(TT) - U)/du
            print*, this%Y, U
            error stop "Compute T From U over 50 iterations."
         endif
      end do

   end function ComputeTemperatureFromU

   function ComputeDensity(this) result(density)  !! kg/m^3
      class(type_chem_phases), intent(inout) :: this
      real(dp) :: density

      density = this%P*this%mean_molecular_weight/(R0*this%T)

   end function ComputeDensity

   function ComputePressure(this) result(pressure)  !! pa
      class(type_chem_phases), intent(inout) :: this
      real(dp) :: pressure

      pressure = this%D*this%T*R0/this%mean_molecular_weight

   end function ComputePressure

   function GetDensity(this) result(density)
      class(type_chem_phases), intent(in) :: this
      real(dp) :: density

      density = this%D

   end function GetDensity

   function GetX(this) result(X)
      class(type_chem_phases), intent(in) :: this
      real(dp), allocatable :: X(:)

      allocate (X(this%n_species))
      X = this%X

   end function GetX

   function GetY(this) result(Y)
      class(type_chem_phases), intent(in) :: this
      real(dp), allocatable :: Y(:)

      allocate (Y(this%n_species))
      Y = this%X*this%species%molar_mass/this%mean_molecular_weight

   end function GetY

   function GetMolarDensity(this) result(molar_density) !! kmol/m3
      class(type_chem_phases), intent(in) :: this
      real(dp), allocatable :: molar_density(:)

      molar_density = this%Y/this%mean_molecular_weight*this%get_density()

   end function GetMolarDensity

   function GetCreationProductionRates(this) result(R)
      class(type_chem_phases), intent(in) :: this
      real(dp), allocatable :: R(:)

      integer :: l
      allocate (R(this%n_species))
      R = 0.0_dp
      do l = 1, this%n_reactions
         R = R + this%reactions(l)%get_creation_production_rates(this%get_molar_density(), this%T, this%Tv)
      end do

   end function GetCreationProductionRates

   function GetDestructionProductionRates(this) result(R)
      class(type_chem_phases), intent(in) :: this
      real(dp), allocatable :: R(:)

      integer :: l
      allocate (R(this%n_species))
      R = 0.0_dp
      do l = 1, this%n_reactions
         R = R + this%reactions(l)%get_desctruction_production_rates(this%get_molar_density(), this%T, this%Tv)
      end do

   end function GetDestructionProductionRates

   function GetNetProductionRates(this) result(R)
      class(type_chem_phases), intent(in) :: this
      real(dp), allocatable :: R(:)

      integer :: l
      allocate (R(this%n_species))
      R = 0.0_dp
      do l = 1, this%n_reactions
         R = R + this%reactions(l)%get_net_production_rates(this%get_molar_density(), this%T, this%Tv)
      end do

   end function GetNetProductionRates

   subroutine GetNetProductionRatesDX(this, R, rr)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(inout) :: R(:)
      real(dp), intent(inout) :: rr(:, :)

      integer :: l
      real(dp), allocatable :: Rt(:)
      real(dp), allocatable :: rrt(:, :)

      allocate(Rt(this%n_species), rrt(this%n_species, this%n_species))
      R = 0.0_dp; rr = 0.0_dp
      Rt = 0.0_dp; rrt = 0.0_dp

      do l = 1, this%n_reactions
         call this%reactions(l)%get_net_production_rates_dx(this%get_molar_density(), this%T, this%Tv, Rt, rrt)
         R = R + Rt
         rr = rr + rrt
      end do

   end subroutine GetNetProductionRatesDX

   function GetFwdRateProgress(this) result(R)
      class(type_chem_phases), intent(in) :: this
      real(dp), allocatable :: R(:)

      integer :: l
      allocate (R(this%n_reactions))
      do l = 1, this%n_reactions
         R(l) = this%reactions(l)%get_fwdrate_progress(this%get_molar_density(), this%T, this%Tv)
      end do

   end function GetFwdRateProgress

   function GetRevRateProgress(this) result(R)
      class(type_chem_phases), intent(in) :: this
      real(dp), allocatable :: R(:)

      integer :: l
      allocate (R(this%n_reactions))
      do l = 1, this%n_reactions
         R(l) = this%reactions(l)%get_revrate_progress(this%get_molar_density(), this%T, this%Tv)
      end do

   end function GetRevRateProgress

   function GetNetRateProgress(this) result(R)
      class(type_chem_phases), intent(in) :: this
      real(dp), allocatable :: R(:)

      integer :: l
      allocate (R(this%n_reactions))
      do l = 1, this%n_reactions
         R(l) = this%reactions(l)%get_netrate_progress(this%get_molar_density(), this%T, this%Tv)
      end do

   end function GetNetRateProgress

   function get_molecular_weights(this) result(molecular_weights) !! kg/kmol
      class(type_chem_phases), intent(in) :: this
      real(dp), allocatable :: molecular_weights(:)

      allocate (molecular_weights(this%n_species))
      molecular_weights = this%species%molar_mass

   end function get_molecular_weights

   function GetMeanMolecularWeight(this) result(mean_molecular_weight) !! kg/kmol
      class(type_chem_phases), intent(in) :: this
      real(dp) :: mean_molecular_weight

      mean_molecular_weight = this%mean_molecular_weight

   end function GetMeanMolecularWeight

   function GetPartialCpR(this, T) result(cp)  !! non-dimensional

      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp), allocatable :: cp(:)

      integer :: i

      allocate (cp(this%n_species))
      do i = 1, this%n_species
         cp(i) = this%species(i)%get_cp_R(optval(T, this%T))
      end do

   end function GetPartialCpR

   function GetPartialCpMole(this, T) result(cp) !! J/(kmol*K)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp), allocatable :: cp(:)

      integer :: i

      allocate (cp(this%n_species))
      do i = 1, this%n_species
         cp(i) = this%species(i)%get_cp_mole(optval(T, this%T))
      end do

   end function GetPartialCpMole

   function GetPartialCpMass(this, T) result(cp) !! J/(kg*K)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp), allocatable :: cp(:)

      integer :: i

      allocate (cp(this%n_species))
      do i = 1, this%n_species
         cp(i) = this%species(i)%get_cp_mass(optval(T, this%T))
      end do
   end function GetPartialCpMass

   function GetCpMole(this, T) result(cp) !! J/(kmol*K)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp) :: cp

      cp = sum(this%X*this%get_partial_cp_mole())
   end function GetCpMole

   function GetCpMass(this, T) result(cp) !! J/(kg*K)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp) :: cp

      cp = sum(this%Y*this%get_partial_cp_mass())

   end function GetCpMass

   function GetCvMole(this, T) result(cv) !! J/(kmol*K)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp) :: cv

      cv = this%get_cp_mole() - R0

   end function GetCvMole

   function GetCvMass(this, T) result(cv) !! J/(kg*K)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp) :: cv

      cv = this%get_cv_mole()/this%get_mean_molecular_weight()

   end function GetCvMass

   function GetPartialEnthalpies_RT(this, T) result(h) !! non-dimensional
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp), allocatable :: h(:)

      integer :: i

      allocate (h(this%n_species))
      do i = 1, this%n_species
         h(i) = this%species(i)%get_enthalpy_RT(optval(T, this%T))
      end do

   end function GetPartialEnthalpies_RT

   function GetPartialEnthalpiesMole(this, T) result(h) !! J/(kmol)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp), allocatable :: h(:)

      integer :: i

      allocate (h(this%n_species))
      do i = 1, this%n_species
         h(i) = this%species(i)%get_enthalpy_mole(optval(T, this%T))
      end do

   end function GetPartialEnthalpiesMole

   function GetPartialEnthalpiesMass(this, T) result(h) !! J/(kg)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp), allocatable :: h(:)

      integer :: i

      allocate (h(this%n_species))
      do i = 1, this%n_species
         h(i) = this%species(i)%get_enthalpy_mass(optval(T, this%T))
      end do

   end function GetPartialEnthalpiesMass

   function GetEnthalpyMole(this, T) result(h)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp) :: h

      h = sum(this%X*this%get_partial_enthalpies_mole(T))

   end function GetEnthalpyMole

   function GetEnthalpyMass(this, T) result(h)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp) :: h

      h = sum(this%Y*this%get_partial_enthalpies_mass(T))

   end function GetEnthalpyMass

   function GetPartialInterEnergies_RT(this, T) result(u) !! non-dimensional
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp), allocatable :: u(:)

      integer :: i

      allocate (u(this%n_species))
      do i = 1, this%n_species
         u(i) = this%species(i)%get_internal_energy_RT(optval(T, this%T))
      end do

   end function GetPartialInterEnergies_RT

   function GetPartialInterEnergiesMole(this, T) result(u) !! J/(kmol)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp), allocatable :: u(:)

      integer :: i

      allocate (u(this%n_species))
      do i = 1, this%n_species
         u(i) = this%species(i)%get_internal_energy_mole(optval(T, this%T))
      end do

   end function GetPartialInterEnergiesMole

   function GetPartialInterEnergiesMass(this, T) result(u) !! J/(kg)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp), allocatable :: u(:)

      integer :: i

      allocate (u(this%n_species))
      do i = 1, this%n_species
         u(i) = this%species(i)%get_internal_energy_mass(optval(T, this%T))
      end do

   end function GetPartialInterEnergiesMass

   function GetInterEnergyMole(this, T) result(u)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp) :: u

      u = sum(this%X*this%get_partial_internal_energies_mole(T))

   end function GetInterEnergyMole

   function GetInterEnergyMass(this, T) result(u)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp) :: u

      u = sum(this%Y*this%get_partial_internal_energies_mass(T))

   end function GetInterEnergyMass

   function GetPartialEntropies_R(this, T) result(s) !! non-dimensional
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp), allocatable :: s(:)

      integer :: i

      allocate (s(this%n_species))
      do i = 1, this%n_species
         s(i) = this%species(i)%get_entropy_R(optval(T, this%T))
      end do

   end function GetPartialEntropies_R

   function GetPartialGibbs_RT(this, T) result(g) !! non-dimensional
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp), allocatable :: g(:)

      integer :: i

      allocate (g(this%n_species))
      do i = 1, this%n_species
         g(i) = this%species(i)%get_gibbs_RT(optval(T, this%T))
      end do

   end function GetPartialGibbs_RT

   function GetGamma(this) result(gamma)
      class(type_chem_phases), intent(in) :: this
      real(dp) :: gamma

      ! gamma = (this%get_enthalpy_mole()-sum(this%X*this%species%thermo%h0))/ &
      !         (this%get_internal_energy_mole()-sum(this%X*this%species%thermo%h0))
      ! gamma = (this%get_enthalpy_mole())/ &
      !         (this%get_internal_energy_mole())
      gamma = (this%get_cp_mole())/ (this%get_cv_mole())

   end function GetGamma

   function GetPartialViscosity(this, T) result(mu) !! Pa-s
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp), allocatable :: mu(:)

      integer :: i

      allocate (mu(this%n_species))
      do i = 1, this%n_species
         mu(i) = this%species(i)%get_viscosity(optval(T, this%T))
      end do

   end function GetPartialViscosity

   function GetPartialThermalConductivity(this, T) result(kappa) !! W/(m K)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp), allocatable :: kappa(:)

      integer :: i

      allocate (kappa(this%n_species))
      do i = 1, this%n_species
         kappa(i) = this%species(i)%get_thermal_conductivity(optval(T, this%T))
      end do

   end function GetPartialThermalConductivity

   function GetViscosity(this, T) result(mu) !! Pa-s
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp) :: mu

      mu = this%gas_transport%get_viscosity(this%X, this%get_molecular_weights(), this%get_partial_viscosity(optval(T, this%T)))

   end function GetViscosity

   function GetThermalConductivity(this, T) result(kappa) !! W/(m K)
      class(type_chem_phases), intent(in) :: this
      real(dp), intent(in), optional :: T
      real(dp) :: kappa

      kappa = this%gas_transport%get_thermal_conductivity( &
              this%X, this%get_molecular_weights(), this%get_partial_thermal_conductivity(optval(T, this%T)))

   end function GetThermalConductivity

   function GetSoundSpeed(this) result(c)
      class(type_chem_phases), intent(in) :: this
      real(dp) :: c

      c = sqrt(this%get_gamma()*this%get_pressure()/this%get_density())

   end function GetSoundSpeed

end module mod_chem_phases
